from utopia.monte_carlo.montecarlo_analysis_paper import run_mc_analysis_paper, set_emission, sim_to_dataframe, reshape_mc_results
from utopia.monte_carlo.monte_carlo_helpers import *
from utopia.utopia import utopiaModel
import numpy as np
import scipy.stats as st
import pandas as pd

data_data = utopiaModel.load_json_file("data/default_data.json")
config_data = utopiaModel.load_json_file("data/default_config.json")

import numpy as np
# To reproduce Rakesh results, use these following discrete values to then generate distributions (lognormal and uniform) from which to sample:

# Plastic density (mpdensity): [900, 1100, 1500] kg/m³ 
# Fragmentation style (FI): [0, 0.3, 0.5, 0.7, 1.0] - ranging from more erosive to more explosive fragmentation 
# Discorporation half time (thalf_deg): [316,1000,10000, 31600, 100000] days 
# Fragmentation timescale (thalf_frag): [10, 31.6, 100, 316] days
# your discrete support values

vals_density = np.array([900, 1100, 1500])
vals_FI = np.array([0, 0.3, 0.5, 0.7, 1.0])
vals_deg = np.array([316, 1000, 10000, 31600, 100000])
vals_frag = np.array([10, 31.6, 100, 316])

def lognorm_params_from_values(values, sigma):
    """
    Return (s, scale) for scipy.stats.lognorm given target values.
    sigma = shape parameter (spread in log space).
    We center the lognormal at the geometric mean of the values.
    """
    mu = np.mean(np.log(values))  # geometric mean in log-space
    scale = np.exp(mu - (sigma**2) / 2)
    return sigma, scale

# choose sigma to control spread (tune manually!)
s_deg, scale_deg = lognorm_params_from_values(vals_deg, sigma=0.5)
s_frag, scale_frag = lognorm_params_from_values(vals_frag, sigma=0.8)
s_density, scale_density = lognorm_params_from_values(vals_density, sigma=0.1)
s_FI, scale_FI = lognorm_params_from_values(vals_FI[vals_FI > 0], sigma=0.5)  # avoid zero for lognorm

param_distributions = {
    "MPdensity_kg_m3": ("lognorm", {"s": s_density, "scale": scale_density}),
    "FI": ("uniform", {"loc": 0.0, "scale": 1.0}),
    "t_half_deg_free": ("lognorm", {"s": s_deg, "scale": scale_deg}),
    "t_frag_gen_FreeSurfaceWater": ("lognorm", {"s": s_frag, "scale": scale_frag})
}

# Visualise the distributions by plotting

for psrsm in param_distributions:
    plot_distribution_from_dict(psrsm, param_distributions)

# Recent studies have challenged the density of biofilm used here (1388 kg m-3) as being too high. see e.g. # https://doi.org/10.1186/s43591-024-00098-2 
# Consider testing chnages in biofouling density

# Set emission scenarios to run that are compatible for comparison with observations data from Kooi et al.

##############################################
# DEFINE ARTICLE -> EMISSION SCENARIOS
##############################################

overlay_emission_scenarios = {
    "Sea Sediment | Bergmann":["Ocean_Surface_Water","Coast_Surface_Water","Coast_Column_Water"],
    "Sea Surface all | Cai":["Surface_Freshwater","Coast_Surface_Water","Impacted_Soil_Surface"],
    "Sea Surface Fragments | Enders":["Ocean_Surface_Water","Coast_Surface_Water"],
    "Sea Surface Fibers | Enders":["Ocean_Surface_Water","Coast_Surface_Water"],
    "Fragments in water | Eo":["Impacted_Soil_Surface","Surface_Freshwater"],
    "Fibers in water | Eo":["Impacted_Soil_Surface","Surface_Freshwater"],
    "Fragments in sediment | Eo":["Impacted_Soil_Surface","Surface_Freshwater"],
    "Fibers in sediment | Eo":["Impacted_Soil_Surface","Surface_Freshwater"],
    "Sea Surface all | Erni-Cassola":["Coast_Surface_Water","Impacted_Soil_Surface"],
    "Fresh water sediments | Imhof":["Surface_Freshwater","Impacted_Soil_Surface"],
    "Sea Surface all | Isobe L1":["Surface_Freshwater","Coast_Surface_Water"],
    "Sea Surface all | Isobe L2":["Surface_Freshwater","Coast_Surface_Water"],
    "Sea Surface all | Isobe L3":["Surface_Freshwater","Coast_Surface_Water"],
    "Sea Surface all | Isobe L4":["Surface_Freshwater","Coast_Surface_Water"],
    "FP all | Scheurer":["Impacted_Soil_Surface"],
    "Sea Surface Fragments | Song":["Impacted_Soil_Surface","Surface_Freshwater"],
    "Sea Surface Fibers | Song":["Impacted_Soil_Surface","Surface_Freshwater"],
    "Sea Surface all | Zhang":["Coast_Surface_Water","Surface_Freshwater"]}
    
    
    
    
# DEFINE ARTICLE->COMPARTMENTS MAPPING
#############################################
overlay_compartments_mapping = {
    "Sea Sediment | Bergmann": ["Sediment_Coast", "Sediment_Ocean"],
    "Sea Surface all | Cai": ["Coast_Surface_Water", "Ocean_Surface_Water"],
    "Sea Surface Fragments | Enders": ["Coast_Surface_Water", "Ocean_Surface_Water"],
    "Sea Surface Fibers | Enders": ["Coast_Surface_Water", "Ocean_Surface_Water"],
    "Fragments in water | Eo": ["Bulk_Freshwater", "Surface_Freshwater"],
    "Fibers in water | Eo": ["Bulk_Freshwater", "Surface_Freshwater"],
    "Fragments in sediment | Eo": ["Sediment_Ocean", "Sediment_Coast", "Sediment_Freshwater"],
    "Fibers in sediment | Eo": ["Sediment_Ocean", "Sediment_Coast", "Sediment_Freshwater"],
    "Sea Surface all | Erni-Cassola": ["Coast_Surface_Water", "Ocean_Surface_Water"],
    "Fresh water sediments | Imhof": ["Bulk_Freshwater", "Surface_Freshwater"],
    "Sea Surface all | Isobe L1": ["Coast_Surface_Water", "Ocean_Surface_Water"], 
    "Sea Surface all | Isobe L2": ["Coast_Surface_Water", "Ocean_Surface_Water"], 
    "Sea Surface all | Isobe L3": ["Coast_Surface_Water", "Ocean_Surface_Water"],
    "Sea Surface all | Isobe L4": ["Coast_Surface_Water", "Ocean_Surface_Water"],
    "FP all | Scheurer": ["Impacted_Soil_Surface", "Sediment_Freshwater"], 
    "Sea Surface Fragments | Song": ["Coast_Surface_Water", "Ocean_Surface_Water"],
    "Sea Surface Fibers | Song": ["Coast_Surface_Water", "Ocean_Surface_Water"],
    "Sea Surface all | Zhang": ["Coast_Surface_Water", "Ocean_Surface_Water"]
}

# NOTE!I have ch anged the mapping for the FP all | Scheurer study compared to the original data from Rakesh. Samples were collected in floodplain soils, so the impacted_soil_surface compartment and freshwater sediment compartment have been included in the compartment mapping.might be more appropriate. Since the model does not simulate floodplain soils as a separate compartment.?

# NOTE 2! Changed entry in overlay compatrments mapping for the "Fragments in sediment | Eo" and "Fibers in sediment | Eo" studies to exclude the beaches compartments.

# Set emission scenarios according to the mapping defined above to the run the model for each parameter combination using the Monte Carlo framework and perform the comparison with the observed data according to mapping.

import seaborn as sns
import matplotlib.pyplot as plt
from utopia.monte_carlo.montecarlo_analysis_paper import run_and_plot_top_results

MC_simulations_results_all={}
for art in overlay_emission_scenarios.keys():
    print(f"Running Monte Carlo analysis for article: {art}")
    MC_simulations_results=[]
    for comp in overlay_emission_scenarios[art]:
        # Set emissions to that compartment and include it in base_data file
        data_data["emiss_dict_g_s"]= set_emission(data_data["emiss_dict_g_s"].copy(), comp, "e", 100)
        
        
        #Run Monte Carlo analysis for each emission scenario and targer compartment
        for comp_T in overlay_compartments_mapping[art]:
            emission_comp=comp
            target_comp=comp_T
            sim= run_mc_analysis_paper(
            base_config=config_data,
            base_data=data_data,
            param_distributions=param_distributions,
            n_cases=20,art=art,emission_comp=emission_comp,target_comp=target_comp
        )
        
            # Convert sim results to dataframe
            df_MC_results = sim_to_dataframe(sim)
            
            
            MC_simulations_results.append(df_MC_results)
    MC_simulations_results = pd.concat(MC_simulations_results, ignore_index=True)
    input_params=["MPdensity_kg_m3", "FI", "t_half_deg_free", "t_frag_gen_FreeSurfaceWater"]
    tidy_mc_results=reshape_mc_results(MC_simulations_results, input_params, art)
    results_sorted_art = tidy_mc_results.sort_values("RMSE", ascending=True)
    # Keep e.g. top 20 best-fitting combinations
    top_results_art = results_sorted_art.head(20)
    run_and_plot_top_results(top_results_art[0:2],data_data, config_data)
    

    # Correlation heatmap for the top results
    corr = top_results_art[0:20][["MPdensity_kg_m3", "FI", "t_half_deg_free", "t_frag_gen_FreeSurfaceWater"]].corr()

    plt.figure(figsize=(8,6))
    sns.heatmap(corr, annot=True, cmap="coolwarm", center=0)
    plt.title("Correlation matrix (Top 20 best fits)")
    plt.show()
    
    
    MC_simulations_results_all[art]=top_results_art

Running Monte Carlo analysis for article: Sea Sediment | Bergmann
Running 'UTOPIA_MC_simulation' Monte Carlo simulation with 20/20 cases...
Drawing random samples for 4 input variables via the 'sobol_random' method... Done
Generating cases... Done

Preprocessing cases: 100%|██████████| 20/20 [00:00<?, ? cases/s]
Running cases:  30%|███       | 6/20 [00:14<00:34,  2.48s/ cases]c:\Users\PradoDomercq\AppData\Local\miniconda3\envs\utopia_envP11\Lib\site-packages\pandas\core\arraylike.py:402: RuntimeWarning: divide by zero encountered in log10
  result = getattr(ufunc, method)(*inputs, **kwargs)
Running cases:  45%|████▌     | 9/20 [00:21<00:25,  2.36s/ cases]c:\Users\PradoDomercq\AppData\Local\miniconda3\envs\utopia_envP11\Lib\site-packages\pandas\core\arraylike.py:402: RuntimeWarning: divide by zero encountered in log10
  result = getattr(ufunc, method)(*inputs, **kwargs)
Running cases: 100%|██████████| 20/20 [00:45<00:00,  2.29s/ cases]
Postprocessing cases: 100%|██████████| 20/20 [00:00<00:00, 5997.86 cases/s]

Simulation complete! Runtime: 0:00:45.844982
Running 'UTOPIA_MC_simulation' Monte Carlo simulation with 20/20 cases...
Drawing random samples for 4 input variables via the 'sobol_random' method... Done
Generating cases... Done

Preprocessing cases: 100%|██████████| 20/20 [00:00<00:00, 6665.56 cases/s] 
Running cases:  30%|███       | 6/20 [00:13<00:32,  2.29s/ cases]c:\Users\PradoDomercq\AppData\Local\miniconda3\envs\utopia_envP11\Lib\site-packages\pandas\core\arraylike.py:402: RuntimeWarning: divide by zero encountered in log10
  result = getattr(ufunc, method)(*inputs, **kwargs)
Running cases:  45%|████▌     | 9/20 [00:20<00:26,  2.37s/ cases]c:\Users\PradoDomercq\AppData\Local\miniconda3\envs\utopia_envP11\Lib\site-packages\pandas\core\arraylike.py:402: RuntimeWarning: divide by zero encountered in log10
  result = getattr(ufunc, method)(*inputs, **kwargs)
Running cases: 100%|██████████| 20/20 [00:45<00:00,  2.26s/ cases]
Postprocessing cases: 100%|██████████| 20/20 [00:00<00:00, 30470.79 cases/s]

Simulation complete! Runtime: 0:00:45.169583
Running 'UTOPIA_MC_simulation' Monte Carlo simulation with 20/20 cases...
Drawing random samples for 4 input variables via the 'sobol_random' method... Done
Generating cases... Done

Preprocessing cases: 100%|██████████| 20/20 [00:00<00:00, 2414.06 cases/s]
Running cases:  30%|███       | 6/20 [00:12<00:30,  2.17s/ cases]c:\Users\PradoDomercq\AppData\Local\miniconda3\envs\utopia_envP11\Lib\site-packages\pandas\core\arraylike.py:402: RuntimeWarning: divide by zero encountered in log10
  result = getattr(ufunc, method)(*inputs, **kwargs)
Running cases:  45%|████▌     | 9/20 [00:19<00:23,  2.14s/ cases]c:\Users\PradoDomercq\AppData\Local\miniconda3\envs\utopia_envP11\Lib\site-packages\pandas\core\arraylike.py:402: RuntimeWarning: divide by zero encountered in log10
  result = getattr(ufunc, method)(*inputs, **kwargs)
Running cases: 100%|██████████| 20/20 [00:45<00:00,  2.29s/ cases]
Postprocessing cases: 100%|██████████| 20/20 [00:00<?, ? cases/s]

Simulation complete! Runtime: 0:00:45.858922
Running 'UTOPIA_MC_simulation' Monte Carlo simulation with 20/20 cases...
Drawing random samples for 4 input variables via the 'sobol_random' method... Done
Generating cases... Done

Preprocessing cases: 100%|██████████| 20/20 [00:00<?, ? cases/s]
Running cases:  30%|███       | 6/20 [00:13<00:30,  2.21s/ cases]c:\Users\PradoDomercq\AppData\Local\miniconda3\envs\utopia_envP11\Lib\site-packages\pandas\core\arraylike.py:402: RuntimeWarning: divide by zero encountered in log10
  result = getattr(ufunc, method)(*inputs, **kwargs)
Running cases:  45%|████▌     | 9/20 [00:19<00:23,  2.14s/ cases]c:\Users\PradoDomercq\AppData\Local\miniconda3\envs\utopia_envP11\Lib\site-packages\pandas\core\arraylike.py:402: RuntimeWarning: divide by zero encountered in log10
  result = getattr(ufunc, method)(*inputs, **kwargs)
Running cases: 100%|██████████| 20/20 [00:43<00:00,  2.16s/ cases]
Postprocessing cases: 100%|██████████| 20/20 [00:00<?, ? cases/s]

Simulation complete! Runtime: 0:00:43.172311
Running 'UTOPIA_MC_simulation' Monte Carlo simulation with 20/20 cases...
Drawing random samples for 4 input variables via the 'sobol_random' method... Done
Generating cases... Done

Preprocessing cases: 100%|██████████| 20/20 [00:00<00:00, 6670.86 cases/s]
Running cases:  30%|███       | 6/20 [00:12<00:28,  2.06s/ cases]c:\Users\PradoDomercq\AppData\Local\miniconda3\envs\utopia_envP11\Lib\site-packages\pandas\core\arraylike.py:402: RuntimeWarning: divide by zero encountered in log10
  result = getattr(ufunc, method)(*inputs, **kwargs)
Running cases:  45%|████▌     | 9/20 [00:18<00:22,  2.06s/ cases]c:\Users\PradoDomercq\AppData\Local\miniconda3\envs\utopia_envP11\Lib\site-packages\pandas\core\arraylike.py:402: RuntimeWarning: divide by zero encountered in log10
  result = getattr(ufunc, method)(*inputs, **kwargs)
Running cases: 100%|██████████| 20/20 [00:40<00:00,  2.01s/ cases]
Postprocessing cases: 100%|██████████| 20/20 [00:00<00:00, 7577.10 cases/s]

combined_df = pd.concat(MC_simulations_results_all.values(), ignore_index=True)

combined_df

# Save results to csv
combined_df.to_csv("UTOPIA_paper_results.csv", index=False)

# Read data from CSV if needed
combined_results = pd.read_csv("UTOPIA_paper_results.csv")
df=combined_results.copy()

combined_results_shorted=df.sort_values("RMSE", ascending=True)
combined_results_shorted

from utopia.monte_carlo.montecarlo_analysis_paper import run_and_plot_top_results
run_and_plot_top_results(combined_results_shorted[0:10],data_data, config_data)

import seaborn as sns
import matplotlib.pyplot as plt

# Correlation heatmap for the top results
corr = combined_results_shorted[0:19][["MPdensity_kg_m3", "FI", "t_half_deg_free", "t_frag_gen_FreeSurfaceWater"]].corr()

plt.figure(figsize=(8,6))
sns.heatmap(corr, annot=True, cmap="coolwarm", center=0)
plt.title("Correlation matrix (Top 20 best fits)")
plt.show()

# Pair-wise relationships for the top results
sns.pairplot(
    combined_results_shorted[0:19],
    vars=["MPdensity_kg_m3", "FI", "t_half_deg_free", "t_frag_gen_FreeSurfaceWater"],
    hue="Target_Compartment",
    corner=True
)
plt.suptitle("Parameter relationships for best-fitting runs", y=1.02)
plt.show()

# Read data from CSV if needed
combined_results = pd.read_csv("UTOPIA_paper_results.csv")
df=combined_results.copy()

combined_results_shorted=df.sort_values("RMSE", ascending=True)

combined_results_shorted[0:20]

import seaborn as sns
import matplotlib.pyplot as plt

# Correlation heatmap for the top results
corr = combined_results_shorted[0:19][["MPdensity_kg_m3", "FI", "t_half_deg_free", "t_frag_gen_FreeSurfaceWater"]].corr()

plt.figure(figsize=(8,6))
sns.heatmap(corr, annot=True, cmap="coolwarm", center=0)
plt.title("Correlation matrix (Top 20 best fits)")
plt.show()

# Write up conclusions and next steps

# Best fits corresponds to 2 archetypical types of plastic:
# 1. Low density plastic (LDP) provides best fits to the monitoring data when framgentation index is low (more erosive fragmentation) and and disintegration half-life is short. Fragmentation timescale seems to not play a significant role according to the correlation matrix.
# 2. High density plastic (HDP), provides best fits when fragmentation index is high (more explosive fragmentation) and disintegration half-life is long. Fragmentation timescale seems to not play a significant role according to the correlation matrix.

# We will run scenarios for the two archetipycal plastics and compare results to the articles data to plot best fits against observed data.

# SET UP SCENARIOS FOR THE TWO ARCHETYPICAL PLASTICS
# 1. Low density plastic (LDP) provides best fits to the monitoring data when framgentation index is low (more erosive fragmentation) and and disintegration half-life is short. Fragmentation timescale seems to not play a significant role according to the correlation matrix.

# Set up scenario simulation for LDP and HDP
from utopia.utopia import utopiaModel
# Define modifications (follow the structure provided in the default configuration data (data_data))
modifications = {
    "MPdensity_kg_m3": 940,
    "MP_composition": "LDP",
    'FI': 0.2,
    "t_half_deg_free": 2000,
    "t_frag_gen_FreeSurfaceWater": 50
}

LDP_data=data_data.copy()
LDP_data["MPdensity_kg_m3"]= modifications["MPdensity_kg_m3"]
LDP_data["MP_composition"]= modifications["MP_composition"]
LDP_data["FI"]= modifications["FI"]
LDP_data["t_half_deg_free"]= modifications["t_half_deg_free"]
LDP_data["t_frag_gen_FreeSurfaceWater"]= modifications["t_frag_gen_FreeSurfaceWater"]    
LDP_data

{'MPdensity_kg_m3': 940,
 'MP_composition': 'LDP',
 'shape': 'sphere',
 'MP_form': 'freeMP',
 'FI': 0.2,
 't_half_deg_free': 2000,
 'heter_deg_factor': 10,
 'biof_deg_factor': 0.5,
 'factor_deepWater_soilSurface': 10,
 'factor_sediment': 100,
 't_frag_gen_FreeSurfaceWater': 50,
 'biof_frag_factor': 2,
 'heter_frag_factor': 100,
 'emiss_dict_g_s': {'Ocean_Surface_Water': {'a': 0,
   'b': 0,
   'c': 0,
   'd': 0,
   'e': 100},
  'Ocean_Mixed_Water': {'a': 0, 'b': 0, 'c': 0, 'd': 0, 'e': 0},
  'Ocean_Column_Water': {'a': 0, 'b': 0, 'c': 0, 'd': 0, 'e': 0},
  'Coast_Surface_Water': {'a': 0, 'b': 0, 'c': 0, 'd': 0, 'e': 0},
  'Coast_Column_Water': {'a': 0, 'b': 0, 'c': 0, 'd': 0, 'e': 0},
  'Surface_Freshwater': {'a': 0, 'b': 0, 'c': 0, 'd': 0, 'e': 0},
  'Bulk_Freshwater': {'a': 0, 'b': 0, 'c': 0, 'd': 0, 'e': 0},
  'Sediment_Freshwater': {'a': 0, 'b': 0, 'c': 0, 'd': 0, 'e': 0},
  'Sediment_Ocean': {'a': 0, 'b': 0, 'c': 0, 'd': 0, 'e': 0},
  'Sediment_Coast': {'a': 0, 'b': 0, 'c': 0, 'd': 0, 'e': 0},
  'Beaches_Soil_Surface': {'a': 0, 'b': 0, 'c': 0, 'd': 0, 'e': 0},
  'Beaches_Deep_Soil': {'a': 0, 'b': 0, 'c': 0, 'd': 0, 'e': 0},
  'Background_Soil_Surface': {'a': 0, 'b': 0, 'c': 0, 'd': 0, 'e': 0},
  'Background_Soil': {'a': 0, 'b': 0, 'c': 0, 'd': 0, 'e': 0},
  'Impacted_Soil_Surface': {'a': 0, 'b': 0, 'c': 0, 'd': 0, 'e': 0},
  'Impacted_Soil': {'a': 0, 'b': 0, 'c': 0, 'd': 0, 'e': 0},
  'Air': {'a': 0, 'b': 0, 'c': 0, 'd': 0, 'e': 0}}}

import pandas as pd
from utopia.utopia import utopiaModel
from utopia.results_processing.process_results import *
from scipy.stats import spearmanr
from utopia.monte_carlo.montecarlo_analysis_paper import extract_log_rel_abundance
from utopia.monte_carlo.montecarlo_analysis_paper import plot_obs_vs_model_with_metrics_True


# Run the model for each emission scenario and target compartment
LDP_simulations_results_all={}
for art in overlay_emission_scenarios.keys():
    print(f"Running Monte Carlo analysis for article: {art}")
    LDP_simulations_results=[]
    for comp in overlay_emission_scenarios[art]:
        # Set emissions to that compartment and include it in base_data file
        LDP_data["emiss_dict_g_s"]= set_emission(LDP_data["emiss_dict_g_s"].copy(), comp, "e", 100)
        for comp_T in overlay_compartments_mapping[art]:
            emission_comp=comp
            target_comp=comp_T
            model = utopiaModel(config=config_data, data=LDP_data)
            model.run()
            processor = ResultsProcessor(model)  # Custom processor to handle results
            processor.estimate_flows()
            processor.generate_flows_dict()
            processor.process_results()
            processor.extract_results_by_compartment()
            df = processor.results_by_comp
            df2 = processor.Results_extended
            # processor.estimate_exposure_indicators()

            # Process model outputs into log-relative abundances per size fraction and compare with observed data
            observed_file = "../src/utopia/data/observed_data_long.xlsx"

            observed_df = pd.read_excel(observed_file)

            observed_df = observed_df.rename(columns={
            "dataset identifier": "Dataset",
            "Article name": "Article",
            "log‑transformed size": "log_Size",
            "log‑transformed abundance": "log_Abundance",
        }
    )
            obs_subset = observed_df[observed_df["Article"] == art]
            x_obs_all = obs_subset["log_Size"].astype(float).values
            y_obs_all = obs_subset["log_Abundance"].astype(float).values

            if len(x_obs_all) > 2:
                slope_obs, intercept_obs = np.polyfit(x_obs_all, y_obs_all, 1)
            else:
                slope_obs, intercept_obs = np.nan, np.nan

            result = {}

            # comp_relative_abundance_resultsaccording to target compartment set in the input parameters

            # Exclude size bins that are not withing the observed data range (0.5 and 5 um)
            excluded_sizes = [0.5, 5]
            filtered_Results = df2[~df2["Size_Fraction_um"].isin(excluded_sizes)]
            rel_abun = extract_log_rel_abundance(
                filtered_Results,
                compartment=target_comp,
                value_col="number_of_particles",
            )

            if rel_abun.empty:
                # fill with NaNs
                result.update(
                    {
                        "Spearman_r": np.nan,
                        "Spearman_p": np.nan,
                        "Pass_Spearman": False,
                        "RMSE": np.nan,
                        "R_squared": np.nan,
                        "Model_slope": np.nan,
                        "Model_intercept": np.nan,
                        "Target_compartment": target_comp,
                        "Emission_compartment": emission_comp,
                    }
                )

            else:

                x_mod_all = rel_abun["log_size"].values
                y_mod_all = rel_abun["log_rel_abundance"].values

                # Spearman correlation
                r_value, p_value = spearmanr(x_mod_all, y_mod_all)

                spearman_threshold = -1

                # Initialize metrics
                rmse, r2, slope_mod, intercept_mod = np.nan, np.nan, np.nan, np.nan
                if np.isnan(r_value):
                    pass_spearman = False
                elif (r_value) == spearman_threshold:
                    pass_spearman = True
                else:
                    pass_spearman = False

                # Compute RMSE, R², slope, intercept only if threshold passed

                if pass_spearman == True and len(x_mod_all) > 2 and len(x_obs_all) > 2:
                    slope_mod, intercept_mod = np.polyfit(x_mod_all, y_mod_all, 1)

                    # Align with observed range
                    x_min = max(x_obs_all.min(), x_mod_all.min())
                    x_max = min(x_obs_all.max(), x_mod_all.max())
                    mask_obs = (x_obs_all >= x_min) & (x_obs_all <= x_max)
                    x_obs_filt = x_obs_all[mask_obs]
                    y_obs_filt = y_obs_all[mask_obs]

                    if len(x_obs_filt) > 2:
                        pred_y = slope_mod * x_obs_filt + intercept_mod
                        ss_res = np.sum((y_obs_filt - pred_y) ** 2)
                        ss_tot = np.sum((y_obs_filt - np.mean(y_obs_filt)) ** 2)
                        r2 = 1 - ss_res / ss_tot if ss_tot > 0 else np.nan
                        rmse = np.sqrt(ss_res / len(x_obs_filt))

                        # Plot observed vs model with metrics only if spearman test is true
                        plot_obs_vs_model_with_metrics_True(
                            x_obs_all,
                            y_obs_all,
                            x_mod_all,
                            y_mod_all,
                            comp_T,
                            art,
                            rmse,
                            r2,
                            MP_density=model.MPdensity_kg_m3,
                            FI=model.FI,
                            t_half_deg=model.t_half_deg_free,
                            t_frag=model.t_frag_gen_FreeSurfaceWater,
                        )

                        # Store all metrics in result dict
                        result.update(
                            {
                                "Spearman_r": r_value,
                                "Spearman_p": p_value,
                                "Pass_Spearman": pass_spearman,
                                "RMSE": rmse,
                                "R_squared": r2,
                                "Model_slope": slope_mod,
                                "Model_intercept": intercept_mod,
                                "Target_compartment": target_comp,
                                "Emission_compartment": emission_comp,
                            }
                        )
                    else:
                        # If not enough points to compute metrics, store NaNs
                        result.update(
                            {
                                "Spearman_r": np.nan,
                                "Spearman_p": np.nan,
                                "Pass_Spearman": False,
                                "RMSE": np.nan,
                                "R_squared": np.nan,
                                "Model_slope": np.nan,
                                "Model_intercept": np.nan,
                                "Target_compartment": target_comp,
                                "Emission_compartment": emission_comp,
                            }
                        )
                else:
                    # Store all metrics in result dict
                    result.update(
                        {
                            "Spearman_r": np.nan,
                            "Spearman_p": np.nan,
                            "Pass_Spearman": False,
                            "RMSE": np.nan,
                            "R_squared": np.nan,
                            "Model_slope": np.nan,
                            "Model_intercept": np.nan,
                            "Target_compartment": target_comp,
                            "Emission_compartment": emission_comp,
                        }
                    )
            LDP_simulations_results.append(result)
    LDP_simulations_results_all[art]=LDP_simulations_results

Running Monte Carlo analysis for article: Sea Sediment | Bergmann

c:\Users\PradoDomercq\AppData\Local\miniconda3\envs\utopia_envP11\Lib\site-packages\pandas\core\arraylike.py:402: RuntimeWarning: divide by zero encountered in log10
  result = getattr(ufunc, method)(*inputs, **kwargs)
c:\Users\PradoDomercq\AppData\Local\miniconda3\envs\utopia_envP11\Lib\site-packages\pandas\core\arraylike.py:402: RuntimeWarning: divide by zero encountered in log10
  result = getattr(ufunc, method)(*inputs, **kwargs)
c:\Users\PradoDomercq\AppData\Local\miniconda3\envs\utopia_envP11\Lib\site-packages\pandas\core\arraylike.py:402: RuntimeWarning: divide by zero encountered in log10
  result = getattr(ufunc, method)(*inputs, **kwargs)
c:\Users\PradoDomercq\AppData\Local\miniconda3\envs\utopia_envP11\Lib\site-packages\pandas\core\arraylike.py:402: RuntimeWarning: divide by zero encountered in log10
  result = getattr(ufunc, method)(*inputs, **kwargs)
c:\Users\PradoDomercq\AppData\Local\miniconda3\envs\utopia_envP11\Lib\site-packages\pandas\core\arraylike.py:402: RuntimeWarning: divide by zero encountered in log10
  result = getattr(ufunc, method)(*inputs, **kwargs)
c:\Users\PradoDomercq\AppData\Local\miniconda3\envs\utopia_envP11\Lib\site-packages\pandas\core\arraylike.py:402: RuntimeWarning: divide by zero encountered in log10
  result = getattr(ufunc, method)(*inputs, **kwargs)

Running Monte Carlo analysis for article: Sea Surface all | Cai

Running Monte Carlo analysis for article: Sea Surface Fragments | Enders

Running Monte Carlo analysis for article: Sea Surface Fibers | Enders

Running Monte Carlo analysis for article: Fragments in water | Eo

Running Monte Carlo analysis for article: Fibers in water | Eo

Running Monte Carlo analysis for article: Fragments in sediment | Eo

c:\Users\PradoDomercq\AppData\Local\miniconda3\envs\utopia_envP11\Lib\site-packages\pandas\core\arraylike.py:402: RuntimeWarning: divide by zero encountered in log10
  result = getattr(ufunc, method)(*inputs, **kwargs)
c:\Users\PradoDomercq\AppData\Local\miniconda3\envs\utopia_envP11\Lib\site-packages\pandas\core\arraylike.py:402: RuntimeWarning: divide by zero encountered in log10
  result = getattr(ufunc, method)(*inputs, **kwargs)
c:\Users\PradoDomercq\AppData\Local\miniconda3\envs\utopia_envP11\Lib\site-packages\pandas\core\arraylike.py:402: RuntimeWarning: divide by zero encountered in log10
  result = getattr(ufunc, method)(*inputs, **kwargs)
c:\Users\PradoDomercq\AppData\Local\miniconda3\envs\utopia_envP11\Lib\site-packages\pandas\core\arraylike.py:402: RuntimeWarning: divide by zero encountered in log10
  result = getattr(ufunc, method)(*inputs, **kwargs)
c:\Users\PradoDomercq\AppData\Local\miniconda3\envs\utopia_envP11\Lib\site-packages\pandas\core\arraylike.py:402: RuntimeWarning: divide by zero encountered in log10
  result = getattr(ufunc, method)(*inputs, **kwargs)
c:\Users\PradoDomercq\AppData\Local\miniconda3\envs\utopia_envP11\Lib\site-packages\pandas\core\arraylike.py:402: RuntimeWarning: divide by zero encountered in log10
  result = getattr(ufunc, method)(*inputs, **kwargs)

Running Monte Carlo analysis for article: Fibers in sediment | Eo

LDP_simulations_results_all

# Flatten it into a list of rows
rows = []
for article, records in LDP_simulations_results_all.items():
    for record in records:
        record["article"] = article  # Add outer key as a column
        rows.append(record)

# Create DataFrame
df = pd.DataFrame(rows)

# 👇 Drop rows with any NaN values
df_clean = df.dropna()

# Reorder columns (optional)
cols = ["article"] + [c for c in df_clean .columns if c != "article"]
df_clean_LDP  = df_clean [cols]

LDP_correlations_shorted=df_clean_LDP.sort_values("RMSE", ascending=True)
LDP_correlations_shorted.head(20)

# Define modifications (follow the structure provided in the default configuration data (data_data))
modifications_HDP = {
    "MPdensity_kg_m3": 1150,
    "MP_composition": "LDP",
    'FI': 0.8,
    "t_half_deg_free": 4000,
    "t_frag_gen_FreeSurfaceWater": 50
}

HDP_data=data_data.copy()
HDP_data["MPdensity_kg_m3"]= modifications_HDP["MPdensity_kg_m3"]
HDP_data["MP_composition"]= modifications_HDP["MP_composition"]
HDP_data["FI"]= modifications_HDP["FI"]
HDP_data["t_half_deg_free"]= modifications_HDP["t_half_deg_free"]
HDP_data["t_frag_gen_FreeSurfaceWater"]= modifications_HDP["t_frag_gen_FreeSurfaceWater"]    
HDP_data

{'MPdensity_kg_m3': 1150,
 'MP_composition': 'LDP',
 'shape': 'sphere',
 'MP_form': 'freeMP',
 'FI': 0.8,
 't_half_deg_free': 4000,
 'heter_deg_factor': 10,
 'biof_deg_factor': 0.5,
 'factor_deepWater_soilSurface': 10,
 'factor_sediment': 100,
 't_frag_gen_FreeSurfaceWater': 50,
 'biof_frag_factor': 2,
 'heter_frag_factor': 100,
 'emiss_dict_g_s': {'Ocean_Surface_Water': {'a': 0,
   'b': 0,
   'c': 0,
   'd': 0,
   'e': 100},
  'Ocean_Mixed_Water': {'a': 0, 'b': 0, 'c': 0, 'd': 0, 'e': 0},
  'Ocean_Column_Water': {'a': 0, 'b': 0, 'c': 0, 'd': 0, 'e': 0},
  'Coast_Surface_Water': {'a': 0, 'b': 0, 'c': 0, 'd': 0, 'e': 0},
  'Coast_Column_Water': {'a': 0, 'b': 0, 'c': 0, 'd': 0, 'e': 0},
  'Surface_Freshwater': {'a': 0, 'b': 0, 'c': 0, 'd': 0, 'e': 0},
  'Bulk_Freshwater': {'a': 0, 'b': 0, 'c': 0, 'd': 0, 'e': 0},
  'Sediment_Freshwater': {'a': 0, 'b': 0, 'c': 0, 'd': 0, 'e': 0},
  'Sediment_Ocean': {'a': 0, 'b': 0, 'c': 0, 'd': 0, 'e': 0},
  'Sediment_Coast': {'a': 0, 'b': 0, 'c': 0, 'd': 0, 'e': 0},
  'Beaches_Soil_Surface': {'a': 0, 'b': 0, 'c': 0, 'd': 0, 'e': 0},
  'Beaches_Deep_Soil': {'a': 0, 'b': 0, 'c': 0, 'd': 0, 'e': 0},
  'Background_Soil_Surface': {'a': 0, 'b': 0, 'c': 0, 'd': 0, 'e': 0},
  'Background_Soil': {'a': 0, 'b': 0, 'c': 0, 'd': 0, 'e': 0},
  'Impacted_Soil_Surface': {'a': 0, 'b': 0, 'c': 0, 'd': 0, 'e': 0},
  'Impacted_Soil': {'a': 0, 'b': 0, 'c': 0, 'd': 0, 'e': 0},
  'Air': {'a': 0, 'b': 0, 'c': 0, 'd': 0, 'e': 0}}}

# Run the model for each emission scenario and target compartment
HDP_simulations_results_all={}
for art in overlay_emission_scenarios.keys():
    print(f"Running Monte Carlo analysis for article: {art}")
    HDP_simulations_results=[]
    for comp in overlay_emission_scenarios[art]:
        # Set emissions to that compartment and include it in base_data file
        HDP_data["emiss_dict_g_s"]= set_emission(HDP_data["emiss_dict_g_s"].copy(), comp, "e", 100)
        for comp_T in overlay_compartments_mapping[art]:
            emission_comp=comp
            target_comp=comp_T
            model = utopiaModel(config=config_data, data=HDP_data)
            model.run()
            processor = ResultsProcessor(model)  # Custom processor to handle results
            processor.estimate_flows()
            processor.generate_flows_dict()
            processor.process_results()
            processor.extract_results_by_compartment()
            df = processor.results_by_comp
            df2 = processor.Results_extended
            # processor.estimate_exposure_indicators()

            # Process model outputs into log-relative abundances per size fraction and compare with observed data
            observed_file = "../src/utopia/data/observed_data_long.xlsx"

            observed_df = pd.read_excel(observed_file)

            observed_df = observed_df.rename(columns={
            "dataset identifier": "Dataset",
            "Article name": "Article",
            "log‑transformed size": "log_Size",
            "log‑transformed abundance": "log_Abundance",
        }
    )
            obs_subset = observed_df[observed_df["Article"] == art]
            x_obs_all = obs_subset["log_Size"].astype(float).values
            y_obs_all = obs_subset["log_Abundance"].astype(float).values

            if len(x_obs_all) > 2:
                slope_obs, intercept_obs = np.polyfit(x_obs_all, y_obs_all, 1)
            else:
                slope_obs, intercept_obs = np.nan, np.nan

            result = {}

            # comp_relative_abundance_resultsaccording to target compartment set in the input parameters

            # Exclude size bins that are not withing the observed data range (0.5 and 5 um)
            excluded_sizes = [0.5, 5]
            filtered_Results = df2[~df2["Size_Fraction_um"].isin(excluded_sizes)]
            rel_abun = extract_log_rel_abundance(
                filtered_Results,
                compartment=target_comp,
                value_col="number_of_particles",
            )

            if rel_abun.empty:
                # fill with NaNs
                result.update(
                    {
                        "Spearman_r": np.nan,
                        "Spearman_p": np.nan,
                        "Pass_Spearman": False,
                        "RMSE": np.nan,
                        "R_squared": np.nan,
                        "Model_slope": np.nan,
                        "Model_intercept": np.nan,
                        "Target_compartment": target_comp,
                        "Emission_compartment": emission_comp,
                    }
                )

            else:

                x_mod_all = rel_abun["log_size"].values
                y_mod_all = rel_abun["log_rel_abundance"].values

                # Spearman correlation
                r_value, p_value = spearmanr(x_mod_all, y_mod_all)

                spearman_threshold = -1

                # Initialize metrics
                rmse, r2, slope_mod, intercept_mod = np.nan, np.nan, np.nan, np.nan
                if np.isnan(r_value):
                    pass_spearman = False
                elif (r_value) == spearman_threshold:
                    pass_spearman = True
                else:
                    pass_spearman = False

                # Compute RMSE, R², slope, intercept only if threshold passed

                if pass_spearman == True and len(x_mod_all) > 2 and len(x_obs_all) > 2:
                    slope_mod, intercept_mod = np.polyfit(x_mod_all, y_mod_all, 1)

                    # Align with observed range
                    x_min = max(x_obs_all.min(), x_mod_all.min())
                    x_max = min(x_obs_all.max(), x_mod_all.max())
                    mask_obs = (x_obs_all >= x_min) & (x_obs_all <= x_max)
                    x_obs_filt = x_obs_all[mask_obs]
                    y_obs_filt = y_obs_all[mask_obs]

                    if len(x_obs_filt) > 2:
                        pred_y = slope_mod * x_obs_filt + intercept_mod
                        ss_res = np.sum((y_obs_filt - pred_y) ** 2)
                        ss_tot = np.sum((y_obs_filt - np.mean(y_obs_filt)) ** 2)
                        r2 = 1 - ss_res / ss_tot if ss_tot > 0 else np.nan
                        rmse = np.sqrt(ss_res / len(x_obs_filt))

                        # Plot observed vs model with metrics only if spearman test is true
                        plot_obs_vs_model_with_metrics_True(
                            x_obs_all,
                            y_obs_all,
                            x_mod_all,
                            y_mod_all,
                            comp_T,
                            art,
                            rmse,
                            r2,
                            MP_density=model.MPdensity_kg_m3,
                            FI=model.FI,
                            t_half_deg=model.t_half_deg_free,
                            t_frag=model.t_frag_gen_FreeSurfaceWater,
                        )

                        # Store all metrics in result dict
                        result.update(
                            {
                                "Spearman_r": r_value,
                                "Spearman_p": p_value,
                                "Pass_Spearman": pass_spearman,
                                "RMSE": rmse,
                                "R_squared": r2,
                                "Model_slope": slope_mod,
                                "Model_intercept": intercept_mod,
                                "Target_compartment": target_comp,
                                "Emission_compartment": emission_comp,
                            }
                        )
                    else:
                        # If not enough points to compute metrics, store NaNs
                        result.update(
                            {
                                "Spearman_r": np.nan,
                                "Spearman_p": np.nan,
                                "Pass_Spearman": False,
                                "RMSE": np.nan,
                                "R_squared": np.nan,
                                "Model_slope": np.nan,
                                "Model_intercept": np.nan,
                                "Target_compartment": target_comp,
                                "Emission_compartment": emission_comp,
                            }
                        )
                else:
                    # Store all metrics in result dict
                    result.update(
                        {
                            "Spearman_r": np.nan,
                            "Spearman_p": np.nan,
                            "Pass_Spearman": False,
                            "RMSE": np.nan,
                            "R_squared": np.nan,
                            "Model_slope": np.nan,
                            "Model_intercept": np.nan,
                            "Target_compartment": target_comp,
                            "Emission_compartment": emission_comp,
                        }
                    )
            HDP_simulations_results.append(result)
    HDP_simulations_results_all[art]=HDP_simulations_results

Running Monte Carlo analysis for article: Sea Sediment | Bergmann

Running Monte Carlo analysis for article: Sea Surface all | Cai

Running Monte Carlo analysis for article: Sea Surface Fragments | Enders

Running Monte Carlo analysis for article: Sea Surface Fibers | Enders

Running Monte Carlo analysis for article: Fragments in water | Eo

Running Monte Carlo analysis for article: Fibers in water | Eo

Running Monte Carlo analysis for article: Fragments in sediment | Eo

Running Monte Carlo analysis for article: Fibers in sediment | Eo

Running Monte Carlo analysis for article: Sea Surface all | Erni-Cassola

Running Monte Carlo analysis for article: Fresh water sediments | Imhof

# Flatten it into a list of rows
rows = []
for article, records in HDP_simulations_results_all.items():
    for record in records:
        record["article"] = article  # Add outer key as a column
        rows.append(record)

# Create DataFrame
df2 = pd.DataFrame(rows)

# 👇 Drop rows with any NaN values
df2_clean = df2.dropna()

# Reorder columns (optional)
cols = ["article"] + [c for c in df2_clean .columns if c != "article"]
df2_clean_HDP  = df2_clean [cols]

HDP_correlations_shorted=df2_clean_HDP.sort_values("RMSE", ascending=True)
HDP_correlations_shorted.head(20)

	Case	Observed_dataset	Emission_Compartment	Target_Compartment	RMSE	R2	Slope	Spearman_r	Pass_Spearman	MPdensity_kg_m3	FI	t_half_deg_free	t_frag_gen_FreeSurfaceWater
0	108	Sea Sediment \| Bergmann	Coast_Column_Water	Sediment_Ocean	0.306841	-0.542379	-2.655127	-1.0	True	1334.859380	0.269348	2665.785177	31.280483
1	68	Sea Sediment \| Bergmann	Coast_Surface_Water	Sediment_Ocean	0.314419	-0.619506	-2.627525	-1.0	True	1334.859380	0.269348	2665.785177	31.280483
2	100	Sea Sediment \| Bergmann	Coast_Column_Water	Sediment_Ocean	0.319249	-0.669642	-2.453298	-1.0	True	1223.294660	0.054891	6510.335737	18.298327
3	112	Sea Sediment \| Bergmann	Coast_Column_Water	Sediment_Ocean	0.336169	-0.851313	-2.519515	-1.0	True	1170.689441	0.163532	10975.160488	41.971349
4	72	Sea Sediment \| Bergmann	Coast_Surface_Water	Sediment_Ocean	0.342700	-0.923945	-2.492947	-1.0	True	1170.689441	0.163532	10975.160488	41.971349
...	...	...	...	...	...	...	...	...	...	...	...	...	...
355	61	Sea Surface all \| Zhang	Surface_Freshwater	Ocean_Surface_Water	3.210394	-19.030002	-2.794831	-1.0	True	1027.531967	0.740818	4154.920691	47.156730
356	47	Sea Surface all \| Zhang	Surface_Freshwater	Coast_Surface_Water	3.554319	-23.551460	-2.733908	-1.0	True	1302.585757	0.539578	9481.499850	12.343862
357	23	Sea Surface all \| Zhang	Coast_Surface_Water	Ocean_Surface_Water	4.112485	-31.867982	-3.125205	-1.0	True	1142.043586	0.902573	3158.382174	79.694563
358	63	Sea Surface all \| Zhang	Surface_Freshwater	Ocean_Surface_Water	4.395308	-36.544219	-3.319479	-1.0	True	1142.043586	0.902573	3158.382174	79.694563
359	37	Sea Surface all \| Zhang	Coast_Surface_Water	Ocean_Surface_Water	4.651851	-41.054845	-3.557464	-1.0	True	1052.671922	0.811533	3073.415520	14.554329

	Case	Observed_dataset	Emission_Compartment	Target_Compartment	RMSE	R2	Slope	Spearman_r	Pass_Spearman	MPdensity_kg_m3	FI	t_half_deg_free	t_frag_gen_FreeSurfaceWater
220	6	Sea Surface all \| Isobe L2	Surface_Freshwater	Coast_Surface_Water	0.192740	0.540131	-0.926349	-1.0	True	938.199911	0.224775	1827.090958	58.014555
240	6	Sea Surface all \| Isobe L3	Surface_Freshwater	Coast_Surface_Water	0.196303	0.512105	-0.926349	-1.0	True	938.199911	0.224775	1827.090958	58.014555
120	97	Fragments in sediment \| Eo	Surface_Freshwater	Sediment_Coast	0.210465	0.830557	-1.757989	-1.0	True	1052.671922	0.811533	3073.415520	14.554329
121	85	Fragments in sediment \| Eo	Surface_Freshwater	Sediment_Coast	0.219214	0.816177	-1.605927	-1.0	True	1067.147583	0.814583	7649.971822	27.484251
260	26	Sea Surface all \| Isobe L4	Surface_Freshwater	Ocean_Surface_Water	0.219473	0.496743	-1.260351	-1.0	True	938.199911	0.224775	1827.090958	58.014555
...	...	...	...	...	...	...	...	...	...	...	...	...	...
219	77	Sea Surface all \| Isobe L1	Coast_Surface_Water	Ocean_Surface_Water	4.229001	-103.626180	-3.557464	-1.0	True	1052.671922	0.811533	3073.415520	14.554329
238	23	Sea Surface all \| Isobe L2	Surface_Freshwater	Ocean_Surface_Water	4.317504	-229.758353	-3.319479	-1.0	True	1142.043586	0.902573	3158.382174	79.694563
358	63	Sea Surface all \| Zhang	Surface_Freshwater	Ocean_Surface_Water	4.395308	-36.544219	-3.319479	-1.0	True	1142.043586	0.902573	3158.382174	79.694563
239	77	Sea Surface all \| Isobe L2	Coast_Surface_Water	Ocean_Surface_Water	4.567080	-257.207700	-3.557464	-1.0	True	1052.671922	0.811533	3073.415520	14.554329
359	37	Sea Surface all \| Zhang	Coast_Surface_Water	Ocean_Surface_Water	4.651851	-41.054845	-3.557464	-1.0	True	1052.671922	0.811533	3073.415520	14.554329

	Case	Observed_dataset	Emission_Compartment	Target_Compartment	RMSE	R2	Slope	Spearman_r	Pass_Spearman	MPdensity_kg_m3	FI	t_half_deg_free	t_frag_gen_FreeSurfaceWater
220	6	Sea Surface all \| Isobe L2	Surface_Freshwater	Coast_Surface_Water	0.192740	0.540131	-0.926349	-1.0	True	938.199911	0.224775	1827.090958	58.014555
240	6	Sea Surface all \| Isobe L3	Surface_Freshwater	Coast_Surface_Water	0.196303	0.512105	-0.926349	-1.0	True	938.199911	0.224775	1827.090958	58.014555
120	97	Fragments in sediment \| Eo	Surface_Freshwater	Sediment_Coast	0.210465	0.830557	-1.757989	-1.0	True	1052.671922	0.811533	3073.415520	14.554329
121	85	Fragments in sediment \| Eo	Surface_Freshwater	Sediment_Coast	0.219214	0.816177	-1.605927	-1.0	True	1067.147583	0.814583	7649.971822	27.484251
260	26	Sea Surface all \| Isobe L4	Surface_Freshwater	Ocean_Surface_Water	0.219473	0.496743	-1.260351	-1.0	True	938.199911	0.224775	1827.090958	58.014555
200	26	Sea Surface all \| Isobe L1	Surface_Freshwater	Ocean_Surface_Water	0.219608	0.717863	-1.260351	-1.0	True	938.199911	0.224775	1827.090958	58.014555
122	95	Fragments in sediment \| Eo	Surface_Freshwater	Sediment_Coast	0.220993	0.813182	-1.619336	-1.0	True	1266.860208	0.753409	4509.321025	33.637734
241	66	Sea Surface all \| Isobe L3	Coast_Surface_Water	Ocean_Surface_Water	0.225292	0.357365	-0.761238	-1.0	True	938.199911	0.224775	1827.090958	58.014555
123	81	Fragments in sediment \| Eo	Surface_Freshwater	Sediment_Coast	0.226836	0.803171	-1.555892	-1.0	True	1027.531967	0.740818	4154.920691	47.156730
261	6	Sea Surface all \| Isobe L4	Surface_Freshwater	Coast_Surface_Water	0.228836	0.452885	-0.926349	-1.0	True	938.199911	0.224775	1827.090958	58.014555
124	91	Fragments in sediment \| Eo	Surface_Freshwater	Sediment_Coast	0.240624	0.778517	-1.524742	-1.0	True	1193.888212	0.679468	7014.614289	66.504082
20	41	Sea Surface all \| Cai	Coast_Surface_Water	Coast_Surface_Water	0.241146	0.833343	-1.488332	-1.0	True	1027.531967	0.740818	4154.920691	47.156730
320	66	Sea Surface Fibers \| Song	Surface_Freshwater	Ocean_Surface_Water	0.242450	0.341708	-1.260351	-1.0	True	938.199911	0.224775	1827.090958	58.014555
140	83	Fibers in sediment \| Eo	Surface_Freshwater	Sediment_Coast	0.246943	-0.238814	-1.236710	-1.0	True	1142.043586	0.902573	3158.382174	79.694563
221	66	Sea Surface all \| Isobe L2	Coast_Surface_Water	Ocean_Surface_Water	0.248879	0.233221	-0.761238	-1.0	True	938.199911	0.224775	1827.090958	58.014555
242	26	Sea Surface all \| Isobe L3	Surface_Freshwater	Ocean_Surface_Water	0.272301	0.061204	-1.260351	-1.0	True	938.199911	0.224775	1827.090958	58.014555
201	6	Sea Surface all \| Isobe L1	Surface_Freshwater	Coast_Surface_Water	0.273270	0.563134	-0.926349	-1.0	True	938.199911	0.224775	1827.090958	58.014555
321	46	Sea Surface Fibers \| Song	Surface_Freshwater	Coast_Surface_Water	0.274691	0.154986	-0.926349	-1.0	True	938.199911	0.224775	1827.090958	58.014555
80	46	Fragments in water \| Eo	Surface_Freshwater	Bulk_Freshwater	0.280510	0.837037	-1.644561	-1.0	True	938.199911	0.224775	1827.090958	58.014555
180	66	Fresh water sediments \| Imhof	Impacted_Soil_Surface	Surface_Freshwater	0.284774	0.219064	-2.776208	-1.0	True	938.199911	0.224775	1827.090958	58.014555

Evaluation of UTOPIA's performance in predicting MPs size dristributionsusing the monaco Monte Carlo framework¶

(https://monaco.readthedocs.io/en/latest/)¶

Import the necessary libraries¶

Load the UTOPIA's default configuration parameters¶

Define the parameter distributions for the Monte Carlo simulation.¶